//============================================================================
// Name        : CppFourierLawHB.cpp
// Author      : Carlos Olivares
// Version     : 0.5
// Copyright   : GPL
// Description : Molecular dynamics simulation
//============================================================================

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string>

#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <time.h>

//#include "Gas.h"
//#include "IOParams.h"
//#include "functions.h"

using namespace std;
#define DIM 3

int indx(int i, int j){
  return (i*DIM+j);
}

int indx(int i, int j, int Nrows){
  return (i*Nrows+j);
}

//DrPhiLJOverr
double DrPhiLJOverr(double delta_r2, double sigma2, double eps){
  //TODO: Make ri6, ri12 global variables.
  double ri6  = (sigma2/delta_r2)*(sigma2/delta_r2)*(sigma2/delta_r2);
  double ri12 = ri6*ri6;
  return -24*eps*(2*ri12 - ri6)/delta_r2;
}

double DrPhiLJ(double delta_r, double sigma2, double eps){
  //TODO: Make ri6, ri12 global variables.
  double delta_r2=delta_r*delta_r;
  double ri6  = (sigma2/delta_r2)*(sigma2/delta_r2)*(sigma2/delta_r2);
  double ri12 = ri6*ri6;
  return -24*eps*(2*ri12 - ri6)/delta_r;
}

// --------------------------------------------
double FSFAxis(double delta_r2, double deltaAxis, double sigma2, double eps, double rc, double DrPhiLJOverr_rc){
  if (delta_r2 < rc*rc){
    return -1.0*(-DrPhiLJOverr(delta_r2, sigma2, eps) + DrPhiLJOverr_rc )*deltaAxis;
  }else{
    return 0.0;
  }
}

double PhiLJ(double delta_r2, double sigma2, double eps){
  double ri6  = (sigma2/delta_r2)*(sigma2/delta_r2)*(sigma2/delta_r2);
  double ri12 = ri6*ri6;
  return 4*eps*(ri12 - ri6);
}

double PhiSF(double delta_r2, double sigma2, double eps, double rc, double PhiLJ_rc, double DrPhiLJOverr_rc){
  double rc2 = rc*rc;
  if (delta_r2 < rc2){
    return PhiLJ(delta_r2, sigma2, eps) - PhiLJ_rc - DrPhiLJOverr_rc*(delta_r2 - rc2)/2;
  }else{
    return 0.0;
  }
}

double PhiWall(double x, double sigma2, double eps,double halfXlen, double rmin){
  double delta_xw;
  double delta_xw2;
  if ( x< -halfXlen){
    delta_xw=(x+halfXlen+rmin);
    delta_xw2=delta_xw*delta_xw;
    return (PhiLJ(delta_xw2, sigma2, eps)+eps);
  }else if ((-halfXlen <= x) and (x<=halfXlen)){
    return 0;
  }else if (halfXlen<x){
    delta_xw=(-x+halfXlen+rmin);
    delta_xw2=delta_xw*delta_xw;
    return PhiLJ(delta_xw2, sigma2, eps)+eps;
  }else{
    cout<<"shit**!!!!!"<<endl;
    return 0;
  }
}

double FWall(double x, double sigma2, double eps,double halfXlen, double rmin){
  double delta_xw;
  if ( x< -halfXlen){
    delta_xw=(x+halfXlen+rmin);
    return -1.0*DrPhiLJ(delta_xw,sigma2,eps);
  }else if ((-halfXlen <= x) and (x<=halfXlen)){
    return 0;
  }else if (halfXlen<x){
    delta_xw=(-x+halfXlen+rmin);
    return  1.0*DrPhiLJ(delta_xw, sigma2, eps);
  }else{
    cout<<"shit!!!!!"<<endl;
    return 0;
  }
}

// Random Number generation with Maxwell distribution
double FOverG(double y){
 return y*exp(-0.5*y*y + y);
}

double randMaxwellDist(gsl_rng *r){
  double Y = 0.0;
  double U = 0.0;
  double c = 2.2039;
  double X = 0.0;

  do {
    Y = gsl_ran_exponential(r,1);
    U = gsl_rng_uniform(r);
  } while(U > FOverG(Y)/c);

  X = Y;
  return X;
}

// --------------------------------------------

class IOParams{
  public:
    double gasDensity;
    int N;
    int Ni;
    double xlen;
    double ylen;
    double T1;
    double T2;
    double dm;
    double otherDeltaR;
    
    string outfilename;
    string outfilename2;
    
    double sqrtT1;
    double sqrtT2;
    int ttime;
    int total_time;
    int timeDim;
    int taux;
    int other_time;
    double dt;
    double dt2;  
    
    
    IOParams(int argc, char **argv);
    ~IOParams();
    
};

class Gas{
  public:
    int N;
    double *vec_r;
    double *vec_v;
    double *mass;

    double xlen,ylen;
    double halfXlen, halfYlen;
    double sigma;
    double rmin;
    
    double eps;
    double sigma2;
    double rc;
    double rc2;
    
    double PhiLJ_rc;
    double DrPhiLJOverr_rc;
    
    double *deltaX;
    double *deltaY;
    double *deltaZ;
    double *deltaR2;
    
    double *Fx;
    double *Fy;
    double *Fz;
    
    double *force;
    double *Ekin_n;
    double *halfV;
    
    double kinEner;
    double potEner;
        
    Gas(int const, double, double);
    Gas(const Gas& orig);
    ~Gas();
    void init1DVec_r(); // FIXME: (TEMPORAL)
    void init1DVec_v(double);
    
    void calculateDelta_lk(int, int);
    void calculateDeltas();
    void calculateForce_lk(int, int);
    void calculateForces();
    void calculateTotalForces_i(int);
    void calculateTotalForces();

    void addWallForces();
    void addXLHeatBath(double, double, gsl_rng*);
    void updatePositions(Gas*, IOParams);
    void updateVelocities(Gas*, IOParams );
    void calculateHalfV();
    void calculateEkin_n();
    void calculatePotentialEnergy();
    void calculateKineticEnergy();
    void addWallPotentialEnergy();
    
    void printVec_r();
    void printHalfV();
    void fprint_vec_r(string);
    void fprint_vec_v(string);
    void readfile_vec_r(string);
    void readfile_vec_v(string);
    void printDeltaX();
    void printFx();
    void printForce();
    void alternateMass(double);
    void printMass();
    void copyVelFrom(const Gas*);
    void copyPosFrom(const Gas*);
    void copyMassFrom(const Gas*);
    void copyForceFrom(const Gas*);

};

Gas::Gas(int const N, double xlen, double ylen){
  this->N = N;
  this->vec_r = new double[N*DIM];
  this->vec_v = new double[N*DIM];
  this->xlen = xlen;
  this->ylen = ylen;
  this->halfXlen = 0.5*xlen;
  this->halfYlen = 0.5*ylen;
  this->sigma=1.0;
  this->rmin = pow(2.0, 1.0/6.0)*this->sigma;
  this->mass = new double[N];
  
  this->deltaX = new double[N*N];
  this->deltaY = new double[N*N];
  this->deltaZ = new double[N*N];
  this->deltaR2= new double[N*N];
  
  this->Fx = new double[N*N];
  this->Fy = new double[N*N];
  this->Fz = new double[N*N];
    
  this->force = new double[N*DIM];
  
  this->Ekin_n = new double[N];
  this->halfV  = new double[N*N];
  
  this->kinEner =0.0;
  this->potEner=0.0;
  
  this->eps   = 0.1;
  this->sigma2= this->sigma*this->sigma;
  this->rc    = 2.5;
  this->rc2   = this->rc*this->rc;

  this->PhiLJ_rc   = PhiLJ(this->rc2, this->sigma2, this->eps);
  this->DrPhiLJOverr_rc = DrPhiLJOverr(this->rc2, this->sigma2, this->eps);
  
  // Initiate vectors Row-major array
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      this->vec_r[indx(i,j)] = 0.0;
      this->vec_v[indx(i,j)] = 0.0;
      this->force[indx(i,j)] = 0.0;
    }
  }
  
  for (int l=0; l < this->N; l++){
    for (int k=0; k < this->N; k++){
      this->Fx[indx(l,k,this->N)] = 0.0;
      this->Fy[indx(l,k,this->N)] = 0.0;
      this->Fz[indx(l,k,this->N)] = 0.0;
    }
  }
  for (int i=0; i< this->N; i++){
    this->mass[i] =1.0;
  }
  
  for (int i=0; i< this->N-1;i++){
    for (int j=0; j<DIM; j++){
      this->vec_r[indx(i,j)] = 0.0;
      this->vec_v[indx(i,j)] = 0.0;
    }
  }
  
  int indx2=0;
  for (int i=0; i < this->N; i++){
    for (int j=0; j < this->N; j++){
      indx2 = i*this->N+j;
      this->deltaX[indx2] = 0.0;
      this->deltaY[indx2] = 0.0;
      this->deltaZ[indx2] = 0.0;
      this->deltaR2[indx2]= 0.0;
    }
  }
  
};

Gas::Gas(const Gas& source){
  this->N = source.N;
  this->vec_r = new double[this->N*DIM];
  this->vec_v = new double[this->N*DIM];
  this->xlen = source.xlen;
  this->ylen = source.ylen;
  this->halfXlen = source.halfXlen;
  this->halfYlen = source.halfYlen;
  this->sigma= source.sigma;
  this->rmin = source.sigma;
  this->mass = new double[this->N];
  
  this->deltaX = new double[this->N*this->N];
  this->deltaY = new double[this->N*this->N];
  this->deltaZ = new double[this->N*this->N];
  this->deltaR2= new double[this->N*this->N];
  
  this->Fx = new double[this->N*this->N];
  this->Fy = new double[this->N*this->N];
  this->Fz = new double[this->N*this->N];
    
  this->force = new double[this->N*DIM];
  
  this->Ekin_n = new double[this->N];
  this->halfV  = new double[this->N*this->N];
  
  this->kinEner =source.kinEner;
  this->potEner =source.potEner;
  
  this->eps   = source.eps;
  this->sigma2= source.sigma2;
  this->rc    = source.rc;
  this->rc2   = source.rc2;

  this->PhiLJ_rc   = source.PhiLJ_rc;
  this->DrPhiLJOverr_rc = source.DrPhiLJOverr_rc;
  
  this->copyPosFrom(&source);
  this->copyVelFrom(&source);
  this->copyMassFrom(&source);
  this->copyForceFrom(&source);
}

void Gas::init1DVec_r(){
  
  double delta_r = 3.0;//(this->xlen)/(this->N+1);
  int j=0; // just for x axis, y=z=0.0
  if (delta_r > this->sigma){
    int h=this->N/2;
    if (this->N%2==1){
      this->vec_r[indx(0,j)] = -h*delta_r;
    }else{
      this->vec_r[indx(0,j)] = -h*delta_r + delta_r/2.0;
    }  
    for (int i=0; i< this->N-1;i++){
      this->vec_r[indx(i+1,j)] = this->vec_r[indx(i,j)]+delta_r;
    }
  }else{
    cout<<"Density to high, decrease to avoid overlap"<<endl;
  }
};// FIXME: (TEMPORAL) 

void Gas::init1DVec_v(double Temp){
    
  const gsl_rng_type *T; T = gsl_rng_default;
  gsl_rng *r; r = gsl_rng_alloc(T);
  long unsigned int seed;
  seed = 100;//(long unsigned int) time(NULL); // Get the time of the system as seed
  gsl_rng_set(r,seed);
  
  int j=0; // just for vx , vy=vz=0.0
  double velmean=0.0;
  for (int i=0; i<this->N; i++){
    this->vec_v[indx(i,j)] = sqrt(Temp)*gsl_ran_gaussian(r, 1);
    velmean+= this->vec_v[indx(i,j)];
  }
  velmean=velmean/N;
  for (int i=0; i<this->N; i++){
    this->vec_v[indx(i,j)] = this->vec_v[indx(i,j)]-velmean;
  }
  
  this->calculateEkin_n();
  this->calculateKineticEnergy();
  double TempFactor=sqrt(0.5*this->N*Temp/this->kinEner);
  
  cout<<"kinEner="<<kinEner<<endl;
  cout<<"TempFactor="<<TempFactor<<endl;
  for (int i=0; i<this->N; i++){
    this->vec_v[indx(i,j)] = TempFactor*this->vec_v[indx(i,j)];
  }
  this->calculateEkin_n();
  this->calculateKineticEnergy();
  cout<<"kinEnerScaled="<<kinEner<<endl;
  gsl_rng_free (r);
};

void Gas::calculateDelta_lk(int l, int k){
  int indx_lk  = l*this->N+k;
  int indx_kl = k*this->N+l;

  this->deltaX[indx_lk] = this->vec_r[indx(k,0)] - this->vec_r[indx(l,0)];
  this->deltaY[indx_lk] = this->vec_r[indx(k,1)] - this->vec_r[indx(l,1)];
  this->deltaZ[indx_lk] = this->vec_r[indx(k,2)] - this->vec_r[indx(l,2)];
  this->deltaR2[indx_lk]= this->deltaX[indx_lk]*this->deltaX[indx_lk]
                      + this->deltaY[indx_lk]*this->deltaY[indx_lk]
                      + this->deltaZ[indx_lk]*this->deltaZ[indx_lk];
    
  this->deltaX[indx_kl] = -this->deltaX[indx_lk];
  this->deltaY[indx_kl] = -this->deltaX[indx_lk];
  this->deltaZ[indx_kl] = -this->deltaX[indx_lk];
  this->deltaR2[indx_kl]= this->deltaR2[indx_lk];
}

void Gas::calculateDeltas(){
//  int indx2=0;
  for (int l=0; l < this->N; l++){
    for (int k=0; k < l; k++){
//      indx2 = l*this->N +k;
      this->calculateDelta_lk(l,k);
//      this->deltaX[indx2] = this->vec_r[indx(k,0)] - this->vec_r[indx(l,0)];
//      this->deltaY[indx2] = this->vec_r[indx(k,1)] - this->vec_r[indx(l,1)];
//      this->deltaZ[indx2] = this->vec_r[indx(k,2)] - this->vec_r[indx(l,2)];
//      this->deltaR2[indx2]= this->deltaX[indx2]*this->deltaX[indx2]
//                          + this->deltaY[indx2]*this->deltaY[indx2]
//                          + this->deltaZ[indx2]*this->deltaZ[indx2];
      
//      cout<<l<<k<<" "<< indx2<<" " <<this->vec_r[indx(k,0)]<<" "<< this->deltaX[indx2] <<endl;
    }
  }
}

void Gas::calculateForce_lk(int l, int k){
  
  int indx_lk = indx(l,k,this->N);
  int indx_kl = indx(k,l,this->N);
  
  this->Fx[indx_lk] = FSFAxis(this->deltaR2[indx_lk], this->deltaX[indx_lk], 
          this->sigma2, this->eps, this->rc, this->DrPhiLJOverr_rc);
  this->Fy[indx_lk] = FSFAxis(this->deltaR2[indx_lk], this->deltaY[indx_lk], 
          this->sigma2, this->eps, this->rc, this->DrPhiLJOverr_rc);
  this->Fz[indx_lk] = FSFAxis(this->deltaR2[indx_lk], this->deltaZ[indx_lk], 
          this->sigma2, this->eps, this->rc, this->DrPhiLJOverr_rc);

  this->Fx[indx_kl] = -this->Fx[indx_lk];
  this->Fy[indx_kl] = -this->Fy[indx_lk];
  this->Fz[indx_kl] = -this->Fz[indx_lk];

}

void Gas::calculateForces(){
  for (int l=0; l < this->N; l++){
    for (int k=0; k < l; k++){
      this->calculateForce_lk(l,k);
//      this->Fx[indx(l,k,this->N)] = FSFAxis(this->deltaR2[indx(l,k,this->N)], this->deltaX[indx(l,k,this->N)], 
//              this->sigma2, this->eps, this->rc, this->DrPhiLJOverr_rc);
//      this->Fy[indx(l,k,this->N)] = FSFAxis(this->deltaR2[indx(l,k,this->N)], this->deltaY[indx(l,k,this->N)], 
//              this->sigma2, this->eps, this->rc, this->DrPhiLJOverr_rc);
//      this->Fz[indx(l,k,this->N)] = FSFAxis(this->deltaR2[indx(l,k,this->N)], this->deltaZ[indx(l,k,this->N)], 
//              this->sigma2, this->eps, this->rc, this->DrPhiLJOverr_rc);
//      
//      this->Fx[indx(k,l,this->N)] = -this->Fx[indx(l,k,this->N)];
//      this->Fy[indx(k,l,this->N)] = -this->Fy[indx(l,k,this->N)];
//      this->Fz[indx(k,l,this->N)] = -this->Fz[indx(l,k,this->N)];
    }
  }
  
//  for (int i=0; i< this->N; i++){
//    this->force[indx(i,0)] = 0.0;
//    this->force[indx(i,1)] = 0.0;
//    this->force[indx(i,2)] = 0.0;
//    for (int k=0; k< this->N; k++){
//      this->force[indx(i,0)] +=  this->Fx[indx(i,k,this->N)];
//      this->force[indx(i,1)] +=  this->Fy[indx(i,k,this->N)];
//      this->force[indx(i,2)] +=  this->Fz[indx(i,k,this->N)];
//    }
//  }
  
}

void Gas::calculateTotalForces_i( int i){
  
  this->force[indx(i,0)] = 0.0;
  this->force[indx(i,1)] = 0.0;
  this->force[indx(i,2)] = 0.0;
  for (int k=0; k< this->N; k++){
    this->force[indx(i,0)] +=  this->Fx[indx(i,k,this->N)];
    this->force[indx(i,1)] +=  this->Fy[indx(i,k,this->N)];
    this->force[indx(i,2)] +=  this->Fz[indx(i,k,this->N)];
  }
}

void Gas::calculateTotalForces(){
  for (int i=0; i< this->N; i++){
    this->calculateTotalForces_i(i);
  }
}

void Gas::addWallForces(){
  int j=0;
  for (int i=0; i<this->N; i++){
    this->force[indx(i,j)] +=  FWall(this->vec_r[indx(i,j)], this->sigma2, this->eps, this->halfXlen,  this->rmin);
  }
}

void Gas::addXLHeatBath(double sqrtT1, double sqrtT2, gsl_rng* r){
  int j=0;
  bool wallflag[this->N];
  for (int i=0; i<this->N; i++){
//    wallflag[i]=false;
    if (this->vec_r[indx(i,j)] < -this->halfXlen){
      this->vec_v[indx(i,j)] = sqrtT1*randMaxwellDist(r);
      this->vec_r[indx(i,j)] = -this->halfXlen;
      wallflag[i]=true;
    }else if(this->halfXlen < this->vec_r[indx(i,j)]){
      this->vec_v[indx(i,j)] = -sqrtT2*randMaxwellDist(r);
      this->vec_r[indx(i,j)] = this->halfXlen;
      wallflag[i]=true;
    }else{
      wallflag[i]=false;
    }
  }
  //TODO: 1. Recalculate the deltaX,Y,Z,R2, Fx, Fy, Fz
  for (int l=0; l<this->N; l++){
    if(wallflag[l]){
      for (int k=0; k<this->N; k++){
        if(k!=l){
          this->calculateDelta_lk(l,k);
          this->calculateForce_lk(l,k);
        }
      }
    }
  }

}

void Gas::updatePositions(Gas* othergas, IOParams io){
    for (int i=0; i< othergas->N; i++){
      for (int j=0; j< DIM; j++){
        this->vec_r[indx(i,j)] = othergas->vec_r[indx(i,j)]
                + io.dt*othergas->vec_v[indx(i,j)]
                + 0.5*io.dt2*othergas->force[indx(i,j)]/othergas->mass[i];
      }
    }
}

void Gas::updateVelocities(Gas* othergas, IOParams io){
for (int i=0; i< this->N; i++){
      for (int j=0; j< DIM; j++){
        this->vec_v[indx(i,j)] = othergas->vec_v[indx(i,j)]
                + 0.5*io.dt*(othergas->force[indx(i,j)] + this->force[indx(i,j)])/this->mass[i];
      }
    }
}

void Gas::calculateHalfV(){
  for(int l=0; l<this->N; l++){
    for (int k=0; k < l; k++){
      this->halfV[indx(l,k,this->N)] = 0.5*PhiSF(this->deltaR2[indx(l,k,this->N)], this->sigma2, 
              this->eps, this->rc, this->PhiLJ_rc, this->DrPhiLJOverr_rc);
      this->halfV[indx(k,l,this->N)] = this->halfV[indx(l,k,this->N)];
    }
  }
}

void Gas::calculateEkin_n(){
  for(int i=0; i< this->N; i++){
    this->Ekin_n[i]= 0.0;
    for(int j=0; j<DIM; j++){
      this->Ekin_n[i]+=this->vec_v[indx(i,j)]*this->vec_v[indx(i,j)];
    }
    this->Ekin_n[i]=0.5*mass[i]*this->Ekin_n[i];
  }
}

void Gas::calculatePotentialEnergy(){
  this->potEner=0.0;
  for(int l=0; l<this->N; l++){
    for (int k=0; k < this->N; k++){
      this->potEner += this->halfV[indx(l,k,this->N)];
    }
  }  
}

void Gas::calculateKineticEnergy(){
  this->kinEner = 0.0;
  for (int i=0; i< this->N; i++){
    this->kinEner += this->Ekin_n[i];
  }
}

void Gas::addWallPotentialEnergy(){
  int j=0;
  for (int i=0; i<this->N; i++){
    this->potEner +=  PhiWall(this->vec_r[indx(i,j)], this->sigma2, this->eps, this->halfXlen,  this->rmin);
  }  
}

Gas::~Gas(){
  delete[] vec_r;
  delete[] vec_v;
  delete[] deltaX;
  delete[] deltaY;
  delete[] deltaZ;
  delete[] deltaR2;
  delete[] Fx;
  delete[] Fy;
  delete[] Fz;
  delete[] mass;
  delete[] force;
  delete[] Ekin_n;
  delete[] halfV;
};

void Gas::printVec_r(){
  
  cout<<"__________"<<endl;
  cout<<"vec_r"<<endl;
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      cout<< this->vec_r[indx(i,j)]<<" ";
    }
    cout<<endl;
  }
  cout<<"__________"<<endl;
};

void Gas::printHalfV(){
  for(int l=0; l<this->N; l++){
    for (int k=0; k < l; k++){
      cout<<this->halfV[indx(l,k,this->N)]<<" ";
    }
    cout<<endl;
  }
}

void Gas::fprint_vec_r(string outfilename){

  ofstream ofile;
  ofile.open(outfilename.c_str());
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      ofile << this->vec_r[indx(i,j)]<<" ";
    }
    ofile << endl;
  }
  ofile.close();
};

void Gas::fprint_vec_v(string outfilename){

  ofstream ofile;
  ofile.open(outfilename.c_str());
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      ofile << this->vec_v[indx(i,j)]<<" ";
    }
    ofile << endl;
  }
  ofile.close();
};

void Gas::readfile_vec_r(string inputfilename){
  ifstream ifile;
  string line;
  ifile.open(inputfilename.c_str());
  
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      ifile >> this->vec_r[indx(i,j)];
    }
  }
  ifile.close();
};

void Gas::readfile_vec_v(string inputfilename){
  ifstream ifile;
  string line;
  ifile.open(inputfilename.c_str());
  
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      ifile >> this->vec_v[indx(i,j)];
    }
  }
  ifile.close();
};

void Gas::printDeltaX(){
for (int l=0; l< this->N; l++){
    for (int k=0; k< this->N; k++){
      cout<< this->deltaX[indx(l, k, this->N)]<<" ";
    }
    cout<<endl;
  }
}

void Gas::printFx(){
for (int l=0; l< this->N; l++){
    for (int k=0; k< this->N; k++){
      cout<< this->Fx[indx(l, k, this->N)]<<" ";
    }
    cout<<endl;
  }
}

void Gas::printForce(){
for (int i=0; i< this->N; i++){
    for (int j=0; j< DIM; j++){
      cout<< this->force[indx(i,j)]<<" ";
    }
    cout<<endl;
  }
}

void Gas::alternateMass(double dm){
  for (int i=0; i<this->N; i++){
    if (i%2==0){
      this->mass[i] =1.0;
    }else{
      this->mass[i] =1.0+dm;
    }
  }
}

void Gas::printMass(){
  for (int i=0; i<this->N; i++){
    cout<<this->mass[i];
    cout<<" ";
  }
  cout<<endl;
}

void Gas::copyVelFrom(const Gas* otherGas){
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      this->vec_v[indx(i,j)]= otherGas->vec_v[indx(i,j)];
    }
  }
}

void Gas::copyPosFrom(const Gas* otherGas){
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
//      cout<<otherGas->N<<endl;
      this->vec_r[indx(i,j)]= otherGas->vec_r[indx(i,j)];
    }
  }
}

void Gas::copyMassFrom(const Gas* otherGas){
  for (int i=0; i<this->N; i++){
      this->mass[i]= otherGas->mass[i];
  }
}

void Gas::copyForceFrom(const Gas* otherGas){
  for (int i=0; i<this->N; i++){
    for (int j=0; j<DIM; j++){
      this->force[indx(i,j)]= otherGas->force[indx(i,j)];
    }
  }
}

/**************************************/
/**************************************/

IOParams::IOParams(int argc, char **argv){
  this->gasDensity=3.0;
  this->N = atoi(argv[1]); //40
  this->Ni = 1./N;
  this->xlen=gasDensity*(N+1);
  this->ylen=4;
  this->T1=atof(argv[2]);//4.0;
  this->T2=atof(argv[3]);//10.0;
  this->dm=atof(argv[4]); // 0.2;
  this->otherDeltaR=atof(argv[5]);
  
  string strN(argv[1]);
  string strT1(argv[2]);
  string strT2(argv[3]);
  string strdm(argv[4]);
  string strODR(argv[5]);
  
  this->outfilename ="N_"+strN+"__T1_"+strT1+"__T2_"+strT2+"__dm_"+strdm+"__dr_"+strODR+".dat";
  this->outfilename2="N_"+strN+"__T1_"+strT1+"__T2_"+strT2+"__dm_"+strdm+"__dr_"+strODR+".trv";
  
  this->sqrtT1=sqrt(T1);
  this->sqrtT2=sqrt(T2);
  this->ttime = 0;
  this->total_time = 4000000; //4000000;
  this->timeDim    = 10000; //40000;
  this->taux = this->total_time/this->timeDim;
  this->other_time = 2000000;

  this->dt=0.001;
  this->dt2 = dt*dt;      
}

IOParams::~IOParams(){

}

/**************************************/
/**************************************/
/**************************************/
/**************************************/
/**************************************/
/**************************************/


int main(int argc, char* argv[]) {
  
  IOParams io(argc, argv);
  ofstream outfile;
  ofstream outfile2;
  
  outfile.open((io.outfilename).c_str());
  outfile2.open((io.outfilename2).c_str());
  // FILE header
  outfile  << "#ttime E_{kin} E_{pot} "<<endl;
  outfile2 << "#ttime ";
  for (int i=0; i<io.N;i++){
    outfile2 << "vx_"<< i <<" ";
    outfile2 << "x_"<< i <<" ";
  }
  outfile2 <<endl;
  const gsl_rng_type *T; T = gsl_rng_default;
  gsl_rng *r; r = gsl_rng_alloc(T);
  long unsigned int seed;
  seed = 100;//(long unsigned int) time(NULL); // Get the time of the system as seed
  gsl_rng_set(r,seed);
  
  // Initial conditions
  Gas gas_old(io.N, io.xlen, io.ylen);
  gas_old.init1DVec_r();
  gas_old.init1DVec_v(io.T1);
  gas_old.alternateMass(io.dm);
  Gas gas_new(gas_old);
  
  gas_old.calculateDeltas();  
  gas_old.calculateForces();
  gas_old.calculateTotalForces();

  while (io.ttime < io.total_time){
    //--- 1. Position update. r(t+dt)
    gas_new.updatePositions(&gas_old, io);    
    
    //--- 2. calculate the forces f(t+dt)
    gas_new.calculateDeltas();
    gas_new.calculateForces();
    gas_new.calculateTotalForces();

    //--- 3. calculate the forces v(t+dt)
    gas_new.updateVelocities(&gas_old, io);
    

    //--- 4. Save data in File.
    if (io.ttime % io.taux == 0){
      outfile  << io.dt*io.ttime << " ";
      outfile2 << io.dt*io.ttime << " ";
      // FIXME: (TEMPORAL)
      for (int i=0; i<gas_new.N; i++){
        outfile2 << gas_new.vec_v[indx(i,0)] << " ";
        outfile2 << gas_new.vec_r[indx(i,0)] << " ";
      }
      outfile2<<endl;
      gas_new.calculateEkin_n();
      gas_new.calculateHalfV();
      gas_new.calculateKineticEnergy();
      gas_new.calculatePotentialEnergy();

      outfile << gas_new.kinEner << " ";
      outfile << gas_new.potEner << " ";
      outfile << endl;
    }
    
    //--- 5. Copy gas_new to gas_old
    gas_old.copyPosFrom(&gas_new);
    gas_old.copyVelFrom(&gas_new);
    gas_old.copyForceFrom(&gas_new);
    io.ttime+=1;
  }  
  outfile.close();
  outfile2.close();
  return 0;
  

//  gas_new.copyPosFrom(&gas_old);
  
//  double const gasDensity=3.0;
//  int const N = atoi(argv[1]); //40
//  int const Ni = 1./N;
//  double xlen=gasDensity*(N+1), ylen=4;
//  double T1=atof(argv[2]);//4.0;
//  double T2=atof(argv[3]);//10.0;
//  double dm=atof(argv[4]); // 0.2;
//  double otherDeltaR=atof(argv[5]);
  
////  string strIniFilename_vec_r(argv[6]);;
////  string strIniFilename_vec_v(argv[7]);
//  
//  double *TempProf = new double[N];
//  
//  string strN(argv[1]);
//  string strT1(argv[2]);
//  string strT2(argv[3]);
//  string strdm(argv[4]);
//  string strODR(argv[5]);
//  
//    
//  
//  string outfilename ="N_"+strN+"__T1_"+strT1+"__T2_"+strT2+"__dm_"+strdm+"__dr_"+strODR+".dat";
//  string outfilename2="N_"+strN+"__T1_"+strT1+"__T2_"+strT2+"__dm_"+strdm+"__dr_"+strODR+".trv";
//  cout<<outfilename<<endl;
//  ofstream outfile;
//  ofstream outfile2;
//  
//  outfile.open(outfilename.c_str());
//  outfile2.open(outfilename2.c_str());
//  // FILE header
//  outfile2 << "#ttime ";
//  for (int i=0; i<N;i++){
//    outfile2 << "x_"<< i <<" ";
//  }
//  for (int i=0; i<N;i++){
//    outfile2 << "vx_"<< i <<" ";
//  }
//  outfile2 << endl;
//  
//

//  // Initial conditions
//  Gas gas_old(N, xlen, ylen); Gas gas_new(N, xlen, ylen);
// 
//    // Mass distribution
//  gas_old.alternateMass(dm);
//  gas_new.alternateMass(dm);
//  
//  gas_old.init1DVec_r(); gas_old.init1DVec_v(T1);
//  //gas_new.init1DVec_r(); gas_new.init1DVec_v(T1);
//  gas_old.printVec_r();
//  //FIXME
////  gas_new.copyVelFrom(gas_old);
////  gas_new.copyPosFrom(gas_old);
//
//  for (int i=0; i<N; i++){
//    gas_new.vec_v[indx(i,0)]= gas_old.vec_v[indx(i,0)];
//  }
//  
//  gas_new.printVec_r();
//  
//  gas_new.printMass();
//  
//  gas_old.calculateDeltas();
//  gas_old.calculateForces();
//
//  const gsl_rng_type *T; T = gsl_rng_default;
//  gsl_rng *r; r = gsl_rng_alloc(T);
//  long unsigned int seed;
//  seed = 100;//(long unsigned int) time(NULL); // Get the time of the system as seed
//  gsl_rng_set(r,seed);
//  
//  double dt=0.0001;
//  double dt2 = dt*dt;  
//
//  // TODO Aqui comienza el While!!!
//  int ttime = 0;
//  int total_time = 4000000;
//  int timeDim    = 40000;
//  int taux = total_time/timeDim;
//  int other_time = 2000000;
//  
//  while (ttime < total_time){
//
//    // Position update.
//    for (int i=0; i< gas_old.N; i++){
//      for (int j=0; j< DIM; j++){
//        gas_new.vec_r[indx(i,j)] = gas_old.vec_r[indx(i,j)]+dt*gas_old.vec_v[indx(i,j)]
//            + 0.5*dt2*gas_old.force[indx(i,j)]/gas_old.mass[i];
//      }
//    }
//    
//    gas_new.calculateDeltas();
//    gas_new.calculateForces();
//    
//    // TODO: LJ Wall potential
//    double rmin= pow(2.0, 1.0/6.0)*gas_new.sigma;
////    cout<<"rmin" << rmin <<endl;
//    double delta_h=0.0, R=0.0, R2=0.0;
//    double wallForce=0.0;
//    double wallEnergy=0.0;
//    bool wallFlag=false;
//    for (int i=0; i< gas_old.N; i++){
//      for (int j=0; j<DIM; j++){
//        wallFlag=false;
//        if (j==0){  //1D case.
//          // Left Wall
//          if (gas_new.vec_r[indx(i,j)]<-gas_new.halfXlen){
//            delta_h = -gas_new.halfXlen - gas_new.vec_r[indx(i,j)]; // positive value
//            if (ttime > other_time){
//              gas_new.vec_v[indx(i,j)] = sqrt(T1)*randMaxwellDist(r);
//              gas_new.vec_r[indx(i,j)] = -gas_new.halfXlen + delta_h;
//              wallEnergy=0.0;
//              wallFlag=true;
//              //TODO in contact with the reservoir doesn't calculate the new velocity from verlet algorithm.
//            }else{
//              R = rmin - delta_h;
//              R2 = R*R;
//              wallForce = -DrPhiLJOverr(R2, gas_new.sigma2, gas_new.eps)*R;
////              cout<<"ttime: " <<ttime<<" "<<"delta_h: "<<delta_h<<" ";
////              cout<<"force: " << gas_new.force[indx(i,j)]<<" ";
//              gas_new.force[indx(i,j)] += wallForce; // force in the positive direction
////              cout<<"+wall: " << gas_new.force[indx(i,j)]<<endl;
//              //TODO: Energy contributions
//              wallEnergy= PhiLJ(R2, gas_new.sigma2, gas_new.eps) +gas_new.eps;
//            }
//          }else if (gas_new.vec_r[indx(i,j)] > gas_new.halfXlen){
//            delta_h = gas_new.vec_r[indx(i,j)]-gas_new.halfXlen; // positive value
//            if (ttime>other_time){
//              gas_new.vec_v[indx(i,j)] = -sqrt(T2)*randMaxwellDist(r); //fabs(gsl_ran_gaussian(r,1));
//              gas_new.vec_r[indx(i,j)] = gas_new.halfXlen - delta_h;
//              wallEnergy=0.0;
//              wallFlag=true;
//            }else{
//              R = rmin - delta_h;
//              R2 = R*R;
//              wallForce = DrPhiLJOverr(R2,gas_new.sigma2,gas_new.eps)*R;
//              gas_new.force[indx(i,j)] +=  wallForce; // force in the negative direction
//              wallEnergy= PhiLJ(R2, gas_new.sigma2, gas_new.eps) +gas_new.eps;
//            }
//          }
//          // Velocity update.
////          cout<<gas_old.vec_r[indx(i,j)]<<" "<<gas_old.vec_v[indx(i,j)]<<endl;
//          if (wallFlag==false)
//            gas_new.vec_v[indx(i,j)] = gas_old.vec_v[indx(i,j)]
//                + 0.5*dt*(gas_old.force[indx(i,j)] + gas_new.force[indx(i,j)])/gas_new.mass[i];
//        }else{
//          gas_new.vec_v[indx(i,j)] = 0.0 ;
//        }
//      }
//    }
//
//    // Save data in File.
//    if (ttime % taux == 0){
//      outfile << dt*ttime << " ";
//      outfile2 << dt*ttime << " ";
//      // FIXME: (TEMPORAL)
//      for (int i=0; i<gas_new.N; i++){
//        outfile2 << gas_new.vec_v[indx(i,0)] << " ";
//        outfile2 << gas_new.vec_r[indx(i,0)] << " ";
//      }
//      outfile2<<endl;
//      gas_new.calculateEkin_n();
//      gas_new.calculateHalfV();
//      gas_new.calculateKineticEnergy();
//      gas_new.calculatePotentialEnergy();
//      gas_new.potEner += wallEnergy;
//      outfile << gas_new.kinEner << " ";
//      outfile << gas_new.potEner << " ";
//      outfile << endl;
////      cout<<".";
//    }
//    
//    for (int i=0; i<gas_old.N; i++){
//      for (int j=0; j<DIM; j++){
//        gas_old.vec_r[indx(i,j)] = gas_new.vec_r[indx(i,j)];
//        gas_old.vec_v[indx(i,j)] = gas_new.vec_v[indx(i,j)];
//        gas_old.force[indx(i,j)] = gas_new.force[indx(i,j)];
//      }
//    }
//    ttime++;
//  }
//  cout<<endl;
//  cout<< "#+++++++++++++++++++" <<endl;
//  outfile.close();
//  outfile2.close();
//  gas_old.printVec_r();
//  gsl_rng_free (r);
//  return 0;
}

